﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Linq.Expressions;
using System.Text;
using System.Threading.Tasks;


using MathNet.Numerics.LinearAlgebra;
using MathNet.Numerics.LinearAlgebra.Double;
using MathNet.Numerics.Statistics;
using MathNet.Symbolics;
using Expr = MathNet.Symbolics.SymbolicExpression;

using NumericalAnalysis.DataInterpolation;
//using static System.Runtime.InteropServices.JavaScript.JSType;

namespace NumericalAnalysis.DataInterpolation
{
    /// <summary>
    /// 拉格朗日插值法  mathnet symbolic Numerics
    /// using MathNet.Symbolics;
    /// using Expr = MathNet.Symbolics.SymbolicExpression;
    /// </summary>
    public class LagrangeInterpolation
    {
        static Expr polynomial = 0;   //多项式表达式
        List<double> polyCoefficients = new List<double>(); //各项系数 幂次高至低
        List<int> coefficienOrders = new List<int>(); //各项幂次

        static Vector<double> theVec = Vector<double>.Build.Dense(2);
        List<Vector<double>> vectors = new List<Vector<double>>();

        int count = 0;

        /// <summary>
        /// 初始化
        /// </summary>
        /// <param name="_vectors"></param>
        public LagrangeInterpolation(IEnumerable<Vector<double>> _vectors)
        {
            count = _vectors.Count();
            vectors = _vectors.ToList();
        }

        public LagrangeInterpolation(IEnumerable<Tuple<double, double>> _tuples)
        {
            count = _tuples.Count();
            vectors = _tuples.Select(v =>
            {
                var vec = Vector<double>.Build.Dense(2);
                vec[0] = v.Item1;
                vec[1] = v.Item2;
                return vec;
            }).ToList();
        }

        /// <summary>
        /// 常规公式解法
        /// </summary>
        /// <returns></returns>
        public Expr FitInterpolation()
        {
            //polynomial = 0;
            if (polynomial.ToString().Contains("t")) return polynomial;

            var t = Expr.Variable("t");

            for (int i = 0; i < count; i++)
            {
                Expr basisFunc = vectors[i][1];

                for (int j = 0; j < i; j++)
                {
                    basisFunc *= (t - vectors[j][0]) / (vectors[i][0] - vectors[j][0]);
                }
                for (int j = i + 1; j < count; j++)
                {
                    basisFunc *= (t - vectors[j][0]) / (vectors[i][0] - vectors[j][0]);
                }
                polynomial += basisFunc.Expand();
            }

            //显示多项式
            Console.WriteLine("--------------------------------");
            Console.WriteLine("多项式 polynomial:");
            Console.WriteLine(polynomial);

            //Console.WriteLine("--------------------------------");
            //Console.WriteLine("多项式简化 polynomial-ExponentialSimplify:");
            //Console.WriteLine(polynomial.ExponentialSimplify());

            //Console.WriteLine("--------------------------------");
            //Console.WriteLine("多项式展开式 polynomial-Expand:");
            //Console.WriteLine(polynomial.Expand());

            polynomial = polynomial.Expand();

            return polynomial;
        }

        /// <summary>
        /// 范德蒙矩阵解法
        /// </summary>
        /// <param name="vectors"></param>
        /// <returns></returns>
        public Expr FitVandermonde()
        {
            if (polynomial.ToString().Contains("t")) return polynomial;
            var t = Expr.Variable("t");

            var vecCount = count;
            DenseMatrix vandermondeMatrix = DenseMatrix.CreateIdentity(vecCount); //new DenseMatrix(vecCount);//3维方阵
            Vector<double> yVec = Vector<double>.Build.Dense(vecCount);

            for (int i = 0; i < count; i++)
            {
                var xValue = vectors[i][0];
                var yValue = vectors[i][1];
                yVec[i] = yValue;
                for (int j = 0; j < vecCount; j++)
                {
                    vandermondeMatrix[i, j] = Math.Pow(xValue, j);
                }
            }

            Vector<double> coefVec = vandermondeMatrix.Inverse() * yVec;

            for (int i = 0; i < coefVec.Count; i++)
            {
                Expr coeffient = coefVec[i];
                Expr item = coeffient;
                for (int j = 0; j < i; j++)
                {
                    item *= t;
                }
                polynomial += item;
            }

            //显示多项式
            Console.WriteLine("--------------------------------");
            Console.WriteLine("多项式 polynomial:");
            Console.WriteLine(polynomial);

            Console.WriteLine("--------------------------------");
            Console.WriteLine("多项式展开式 polynomial-Expand:");
            Console.WriteLine(polynomial.Expand());

            return polynomial;
        }

        /// <summary>
        /// 从给出的向量计算
        /// </summary>
        /// <param name="sourceVecs"></param>
        /// <returns></returns>
        public IList<Vector<double>> CalculateBy(IList<Vector<double>> sourceVecs)
        {
            for (int i = 0; i < sourceVecs.Count(); i++)
            {
                double xValue = sourceVecs[i][0];
                var symbols = new Dictionary<string, FloatingPoint> { { "t", xValue } };
                double yValue = polynomial.Evaluate(symbols).RealValue;
                //double   yValue = CalculateBy(xValue);
                sourceVecs[i][1] = yValue;
            }
            return sourceVecs;
        }

        /// <summary>
        /// 从给出的数列计算
        /// </summary>
        /// <param name="xValues"></param>
        /// <returns></returns>
        public IList<double> CalculateBy(IList<double> xValues)
        {
            List<double> yValues = new List<double>();
            for (int i = 0; i < xValues.Count(); i++)
            {
                var symbols = new Dictionary<string, FloatingPoint> { { "t", xValues[i] } };
                double yValue = polynomial.Evaluate(symbols).RealValue;

                yValues.Add(yValue);
            }
            return yValues;
        }

        /// <summary>
        /// 从给出的数计算
        /// </summary>
        /// <param name="xValue"></param>
        /// <returns></returns>
        private double CalculateBy(double xValue)
        {
            //复杂多项式避免使用Compile方法计算  将严重拖慢性能
            //Func<double, double> func = polynomial.Compile("t"); 
            //double yValue = func(xValue);

            var symbols = new Dictionary<string, FloatingPoint> { { "t", xValue } };
            double yValue = polynomial.Evaluate(symbols).RealValue;

            return yValue;
        }


    }

    public class LagrangeInterpolationTest
    {
        string var_t = "t";
         
        /// <summary>
        /// 自定义插值点
        /// </summary>
        /// <returns></returns>
        public (List<Tuple<double, double>>, List<Tuple<double, double>>) ExampleHomeWork()
        {  
            List<double> values = new List<double>()
                { 12, 9, 9, 10, 18, 24, 28, 27, 25, 20, 18, 15, 13 }; 
            List<Tuple<double, double>> interpPnts = new List<Tuple<double, double>>();
            for (int i = 0; i < values.Count; i++)
            {
                double xk = i * 2;
                double yk = values[i];
                Tuple<double, double> tuple = new Tuple<double, double>(xk, yk);
                interpPnts.Add(tuple);
            } 

            List<double> checkXs = NumericalService.LinSpace(0 + 0.1, 24 - 0.1, 0.1);
            ////进行插值   
            var interp = new LagrangeInterpolation(interpPnts);
            //var expression = interp.FitInterpolation();  //插值函数
            var expression = interp.FitVandermonde();  //插值函数
            var checkdPnts = expression.GeneratePointBy(var_t, checkXs);

            return (interpPnts.ToList(), checkdPnts.ToList());
        }


        /// <summary>
        /// 正弦函数
        /// </summary>
        /// <returns></returns>
        public (List<Tuple<double, double>>, List<Tuple<double, double>>) ExampleSin()
        {
            List<double> interpXs = NumericalService.LinSpace(0, Math.PI * 2, Math.PI * 2 * 0.1);
            List<double> checkXs = NumericalService.LinSpace(0, Math.PI * 2, Math.PI * 2 * 0.02);   //构造sin函数插值点 
         
            var t = Expr.Variable(var_t);
            Expr sourceExpr = t.Sin();  //原函数 
            //原函数的离散点
            var interpPnts = sourceExpr.GeneratePointBy(var_t, interpXs).ToList();

            ////进行插值   
            var interp = new LagrangeInterpolation(interpPnts);
            //var expression = interp.FitInterpolation();  //插值函数
            var expression = interp.FitVandermonde();  //插值函数
            var checkdPnts = expression.GeneratePointBy(var_t, checkXs);  

            return (interpPnts.ToList(), checkdPnts.ToList());
        }
         
        /// <summary>
        /// 龙格函数
        /// </summary>
        /// <returns></returns>
        public (List<Tuple<double, double>>, List<Tuple<double, double>>) ExampleRunge()
        {
            List<double> interpXs = NumericalService.LinSpace(-5, 5, 0.5);
            List<double> checkXs = NumericalService.LinSpace(-4, 4, 0.02);
           
            Expr sourceExpr = "1 / (t * t + 1)";
            //原函数的离散点
            var interpPnts = sourceExpr.GeneratePointBy(var_t, interpXs).ToList();

            ////进行插值   
            var interp = new LagrangeInterpolation(interpPnts);
            var expression = interp.FitInterpolation();  //插值函数
            var checkdPnts = expression.GeneratePointBy(var_t, checkXs);

            return (interpPnts.ToList(), checkdPnts.ToList());


        }

    }

}
